# -*- coding: utf-8 -*-
"""
An example python script to calculate urban emissivity

@author: Ciel
"""
import numpy as np
import netCDF4 as nc 

Dir='G:\\SSP5UC_DataPreparation\\deposit\\Surface\\' # The directory you are working with

# Read surface data
SurfData=nc.Dataset(Dir+'b.e21.SSP585UC_IndividualSoil_surfdata_0.9x1.25_16pfts.nc')

# ---------------calculate urban emissivity -------------------
PCT_URBAN = SurfData.variables["PCT_URBAN"][:,:,:]
PCT_URBAN_TOT=np.sum(PCT_URBAN,axis=0) # Sum up the total urban fractions

# Read the variables needed in Eq. S1-S6
WTROAD_PERV = SurfData.variables["WTROAD_PERV"][:,:,:]
WTLUNIT_ROOF = SurfData.variables["WTLUNIT_ROOF"][:,:,:]
EM_WALL = SurfData.variables["EM_WALL"][:,:,:]
EM_ROOF = SurfData.variables["EM_ROOF"][:,:,:]
EM_PERROAD = SurfData.variables["EM_PERROAD"][:,:,:]
EM_IMPROAD = SurfData.variables["EM_IMPROAD"][:,:,:]
CANYON_HWR = SurfData.variables["CANYON_HWR"][:,:,:]

vf_sky_road = np.sqrt(CANYON_HWR**2.0 + 1.0) - CANYON_HWR # the sky-road view factor (Eq.S1)
vf_sky_wall = 0.5 * (CANYON_HWR + 1.0 - np.sqrt(CANYON_HWR**2.0 + 1.0))/CANYON_HWR # the sky-wall view factor (Eq.S2)
em_road = EM_PERROAD* WTROAD_PERV + EM_IMPROAD*(1.0 - WTROAD_PERV) # the emissivity of road (Eq.S3)
em_canyon = em_road*vf_sky_road + EM_WALL*vf_sky_wall # the emissivity of canyon (Eq.S4)
em_urban = EM_ROOF*WTLUNIT_ROOF+ em_canyon*(1.0-WTLUNIT_ROOF) # the urban emissivity (Eq.S5)

# For each grid, we have the urban emissivity for the TBD, HD, and MD urban class (TBD:0;HD:1;MD:2)
# The effective emissivity for the whole urban land unit is the area-weighted mean emissivity of the thre urban classes.
Emi_U=(em_urban[0,:,:]*PCT_URBAN[0,:,:] + em_urban[1,:,:]*PCT_URBAN[1,:,:] + em_urban[2,:,:]*PCT_URBAN[2,:,:])/PCT_URBAN_TOT